Cryptic persistence of male gonad molecular groundplan in parthenogenetic Bacillus species
Abstract:
After the loss of a trait, theory predicts that the molecular machinery underlying its phenotypic expression should decay. Yet, empirical evidence is contrasting. Here, we test the hypotheses that (1) the molecular ground plan of a lost trait could persist due to pleiotropic effects on other traits and (2) that gene co-expression network architecture could constrain individual gene expression. Our testing ground has been the Bacillus stick insect species complex, which contains close relatives that are either bisexual or parthenogenetic. After the identification of genes expressed in male reproductive tissues in a bisexual species, we investigated their gene co-expression network structure in two parthenogenetic species. We found that gene co-expression within the male gonads was preserved in parthenogens. Furthermore, parthenogens did not show relaxed selection on genes upregulated in male gonads in the bisexial species. As these genes were mostly expressed in female gonads, this preservation could be driven by pleiotropic interactions and an ongoing role in female reproduction. Connectivity within the network also played a key role, with highly connected - and more pleiotropic - genes within male gonad also having a gonad-biased expression in parthenogens. Our findings provide novel insight into the mechanisms which could underlie the production of rare males in parthenogenetic lineages; more generally, they provide an example of the cryptic persistence of a lost trait molecular ground plan, driven by gene pleiotropy on other traits and within their co-expression network.
Significance:
Loss of traits commonly occurs in diverse lineages of organisms. Here we investigate what happens to genes and regulatory networks associated with these traits, using parthenogenetic insect species as a model. We investigated the fate of genes and gene regulatory networks associated with male gonads in a bisexual species in closely related parthenogens. Rather than showing signs of disuse and decay, they have been partially preserved in parthenogens. More highly pleiotropic genes in male gonads were more likely to have a gonad-biased expression profile in parthenogens. These results highlight the role of pleiotropy in the cryptic persistence of a trait molecular ground plan, despite its phenotypical absence.
library requirements
require(WGCNA)
require(data.table)
require(tidyverse)
require(grid)
require(gridExtra)
require(ggplot2)
require(dplyr)
require(tidyr)
require(networkD3)
require(ggpubr)
require(gtools)
require(VennDiagram)
require(ggridges)
require(GeneOverlap)
require(ggpubr)
options(stringsAsFactors = FALSE);
import B. grandii expression data
femData = read.csv("BGM_RSEM.TMM.EXPR.matrix", sep =" ", header = TRUE)
dim(femData);## [1] 2785 25names(femData);## [1] "transcript" "BGM_F_G_RNA10" "BGM_F_L_RNA10" "BGM_F_G_RNA12"
## [5] "BGM_F_L_RNA12" "BGM_F_G_RNA14" "BGM_F_L_RNA14" "BGM_F_G_RNA15"
## [9] "BGM_F_L_RNA15" "BGM_F_G_RNA17" "BGM_F_L_RNA17" "BGM_F_G_RNA19"
## [13] "BGM_F_L_RNA19" "BGM_M_G_RNA04" "BGM_M_L_RNA04" "BGM_M_G_RNA06"
## [17] "BGM_M_L_RNA06" "BGM_M_G_RNA08" "BGM_M_L_RNA08" "BGM_M_G_RNA20"
## [21] "BGM_M_L_RNA20" "BGM_M_G_RNA21" "BGM_M_L_RNA21" "BGM_M_G_RNA22"
## [25] "BGM_M_L_RNA22"datExpr0 = as.data.frame(t(femData[, -c(1)]));
names(datExpr0) = femData$transcript;
rownames(datExpr0) = names(femData)[-c(1)];
gsg = goodSamplesGenes(datExpr0, verbose = 3);## Flagging genes and samples with too many missing values...
## ..step 1gsg$allOK## [1] TRUEif (!gsg$allOK)
{
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
# Remove the offending genes and samples from the data:
datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
import B. grandii trait data
traitData = read.csv("traits_description.csv");
traitData = traitData[, c(2,3:8) ];
dim(traitData);## [1] 24 7names(traitData);## [1] "sample" "male_gonads" "female_gonads" "male_legs"
## [5] "female_legs" "sex" "tissue"Samples = rownames(datExpr0);
traitRows = match(Samples, traitData$sample);
datTraits = traitData[traitRows, -1];
rownames(datTraits) = traitData[traitRows, 1];
collectGarbage();
# Plot samples and traits tree
sampleTree2 = hclust(dist(datExpr0), method = "average");
traitColors = numbers2colors(datTraits, signed = TRUE);
plotDendroAndColors(sampleTree2, traitColors,
groupLabels = names(datTraits),
main = "Sample dendrogram and trait heatmap");
infer B. grandii coexpression network
allowWGCNAThreads() ## Allowing multi-threading with up to 4 threads.powers = c(c(1:10), seq(from = 5, to=32, by=2))
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5,networkType = "signed")## pickSoftThreshold: will use block size 2785.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 2785 of 2785
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.2920 3.5200 0.4940 1500.0 1530.0 1730
## 2 2 0.4860 3.7100 0.6660 976.0 930.0 1270
## 3 3 0.4060 1.9200 0.5880 709.0 639.0 1060
## 4 4 0.2240 0.7960 0.3830 551.0 480.0 942
## 5 5 0.0643 0.2870 0.1670 446.0 379.0 855
## 6 5 0.0643 0.2870 0.1670 446.0 379.0 855
## 7 6 0.0128 -0.0923 -0.0322 372.0 307.0 788
## 8 7 0.1680 -0.3040 0.0261 317.0 252.0 733
## 9 7 0.1680 -0.3040 0.0261 317.0 252.0 733
## 10 8 0.3740 -0.4710 0.2080 275.0 210.0 687
## 11 9 0.5080 -0.5650 0.3680 242.0 177.0 647
## 12 9 0.5080 -0.5650 0.3680 242.0 177.0 647
## 13 10 0.5730 -0.6370 0.4550 215.0 152.0 612
## 14 11 0.6580 -0.6810 0.5750 192.0 130.0 582
## 15 13 0.7100 -0.7510 0.6650 158.0 98.3 532
## 16 15 0.7600 -0.7990 0.7550 133.0 75.8 490
## 17 17 0.7900 -0.8350 0.8110 113.0 59.6 455
## 18 19 0.8190 -0.8650 0.8680 98.2 47.2 425
## 19 21 0.8340 -0.8770 0.8960 86.1 38.3 398
## 20 23 0.8370 -0.8950 0.9150 76.3 31.3 376
## 21 25 0.8490 -0.9080 0.9280 68.1 25.8 356
## 22 27 0.8510 -0.9180 0.9370 61.2 21.3 338
## 23 29 0.8700 -0.9220 0.9520 55.4 17.7 322
## 24 31 0.8840 -0.9280 0.9610 50.3 14.7 308sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.80,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
net = blockwiseModules(datExpr0, power = 19,
TOMType = "signed", networkType = "signed", minModuleSize = 30, mergeCutHeight = 0.2, corOptions = list(use = 'p', maxPOutliers = 0.1),
numericLabels = TRUE, pamRespectsDendro = FALSE, replaceMissingAdjacencies = TRUE,
saveTOMs = TRUE, saveTOMFileBase = "test",
verbose = 3, maxBlockSize = 6000, corType = "bicor")## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will use 4 parallel threads.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ..saving TOM for block 1 into file test-block.1.RData
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..removing 18 genes from module 1 because their KME is too low.
## ..removing 3 genes from module 2 because their KME is too low.
## ..removing 5 genes from module 3 because their KME is too low.
## ..removing 4 genes from module 4 because their KME is too low.
## ..removing 1 genes from module 5 because their KME is too low.
## ..removing 5 genes from module 7 because their KME is too low.
## ..removing 3 genes from module 10 because their KME is too low.
## ..reassigning 2 genes from module 6 to modules with higher KME.
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0.2
## Calculating new MEs...# Plot the dendrogram and the module colors underneath
sizeGrWindow(12, 9)
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
# Recalculate MEs with color labels
nGenes = ncol(datExpr0);
nSamples = nrow(datExpr0);
MEs0 = moduleEigengenes(datExpr0, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "all.obs", method = c("pearson"));
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
moduleTraitPvalue <- p.adjust(moduleTraitPvalue, method = "hochberg", n = length(moduleTraitPvalue))
# Display the correlation values within a heatmap plot
sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(90),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.8,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
characterize B. grandii coexpression modules
### module genes
black_module <- names(datExpr0)[moduleColors=="black"]
blue_module <- names(datExpr0)[moduleColors=="blue"]
brown_module <- names(datExpr0)[moduleColors=="brown"]
red_module <- names(datExpr0)[moduleColors=="red"]
green_module <- names(datExpr0)[moduleColors=="green"]
yellow_module <- names(datExpr0)[moduleColors=="yellow"]
turquoise_module <- names(datExpr0)[moduleColors=="turquoise"]
grey_module <- names(datExpr0)[moduleColors=="grey"]
### sizes
black_moduleln <- length(black_module)
blue_moduleln <- length(blue_module)
brown_moduleln <- length(brown_module)
red_moduleln <- length(red_module)
green_moduleln <- length(green_module)
yellow_moduleln <- length(yellow_module)
turquoise_moduleln <- length(turquoise_module)
grey_moduleln <- length(grey_module)
# Define variable containing the variable column of datTrait
mg = as.data.frame(datTraits$male_gonads);
names(mg) = "mg"
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "all.obs", method = c("spearman")));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datExpr0, mg, use = "all.obs", method = c("spearman")));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(mg), sep="");
names(GSPvalue) = paste("p.GS.", names(mg), sep="");
geneModuleMembership_re <- setDT(geneModuleMembership, keep.rownames = TRUE)[]
geneTraitSignificance_re <- setDT(geneTraitSignificance, keep.rownames = TRUE)[]mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_turquoise <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% turquoise_module)
pearson_turquoise <- cor.test(abs(mmembership_gsignificance_turquoise$GS.mg), abs(mmembership_gsignificance_turquoise$MMturquoise), method = "pearson")
turquoise_plot <- ggplot() + geom_point(data=mmembership_gsignificance_turquoise, aes(abs(GS.mg),abs(MMturquoise)), size = 5, alpha = 0.3, shape= 20, color="turquoise") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"turquoise module pval: ", round(pearson_turquoise$p.value,digits = 3), " rho: ", round(pearson_turquoise$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_green <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% green_module)
pearson_green <- cor.test(abs(mmembership_gsignificance_green$GS.mg), abs(mmembership_gsignificance_green$MMgreen), method = "pearson")
green_plot <- ggplot() + geom_point(data=mmembership_gsignificance_green, aes(abs(GS.mg),abs(MMgreen)), size = 5, alpha = 0.3, shape= 20, color="green") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"green module pval: ", round(pearson_green$p.value,digits = 3), " rho: ", round(pearson_green$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_red <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% red_module)
pearson_red <- cor.test(abs(mmembership_gsignificance_red$GS.mg), abs(mmembership_gsignificance_red$MMred), method = "pearson")
red_plot <- ggplot() + geom_point(data=mmembership_gsignificance_red, aes(abs(GS.mg),abs(MMred)), size = 5, alpha = 0.3, shape= 20, color="red") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"red module pval: ", round(pearson_red$p.value,digits = 3), " rho: ", round(pearson_red$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_yellow <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% yellow_module)
pearson_yellow <- cor.test(abs(mmembership_gsignificance_yellow$GS.mg), abs(mmembership_gsignificance_yellow$MMyellow), method = "pearson")
yellow_plot <- ggplot() + geom_point(data=mmembership_gsignificance_yellow, aes(abs(GS.mg),abs(MMyellow)), size = 5, alpha = 0.3, shape= 20, color="yellow") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"yellow module pval: ", round(pearson_yellow$p.value,digits = 3), " rho: ", round(pearson_yellow$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_brown <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% brown_module)
pearson_brown <- cor.test(abs(mmembership_gsignificance_brown$GS.mg), abs(mmembership_gsignificance_brown$MMbrown), method = "pearson")
brown_plot <- ggplot() + geom_point(data=mmembership_gsignificance_brown, aes(abs(GS.mg),abs(MMbrown)), size = 5, alpha = 0.3, shape= 20, color="brown") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"brown module pval: ", round(pearson_brown$p.value,digits = 3), " rho: ", round(pearson_brown$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_black <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% black_module)
pearson_black <- cor.test(abs(mmembership_gsignificance_black$GS.mg), abs(mmembership_gsignificance_black$MMblack), method = "pearson")
black_plot <- ggplot() + geom_point(data=mmembership_gsignificance_black, aes(abs(GS.mg),abs(MMblack)), size = 5, alpha = 0.3, shape= 20, color="black") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"black module pval: ", round(pearson_black$p.value,digits = 3), " rho: ", round(pearson_black$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
mmembership_gsignificance <- merge(geneModuleMembership_re, geneTraitSignificance_re, by="rn")
mmembership_gsignificance_blue <- subset(mmembership_gsignificance, mmembership_gsignificance$rn %in% blue_module)
pearson_blue <- cor.test(abs(mmembership_gsignificance_blue$GS.mg), abs(mmembership_gsignificance_blue$MMblue), method = "pearson")
blue_plot <- ggplot() + geom_point(data=mmembership_gsignificance_blue, aes(abs(GS.mg),abs(MMblue)), size = 5, alpha = 0.3, shape= 20, color="blue") +
labs(x="Gene Significance", y = "module membership", title = paste0(
"blue module pval: ", round(pearson_blue$p.value,digits = 3), " rho: ", round(pearson_blue$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=13,face="bold")) + theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "blue", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
pmodules <- grid.arrange(turquoise_plot,
green_plot,
yellow_plot,
blue_plot,
black_plot,
red_plot,
brown_plot)
pmodules## TableGrob (3 x 3) "arrange": 7 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
## 4 4 (2-2,1-1) arrange gtable[layout]
## 5 5 (2-2,2-2) arrange gtable[layout]
## 6 6 (2-2,3-3) arrange gtable[layout]
## 7 7 (3-3,1-1) arrange gtable[layout]
B. grandii coexpression network plots (Figure 1)
moduleTraitCor_heat <- moduleTraitCor
moduleTraitCor_heat_df <- moduleTraitCor_heat %>% as.data.frame
moduleTraitCor_heat_df <- setDT(moduleTraitCor_heat_df, keep.rownames = TRUE)[]
moduleTraitCor_heat_df %>% gather(trait, value, 2:7) -> moduleTraitCor_heat_df
moduleTraitCor_heat_df["padj"] <- moduleTraitPvalue
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEblack"] <- paste("black - A", black_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEblue"] <- paste("blue - B", blue_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEbrown"] <- paste("brown - C", brown_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEgreen"] <- paste("green - D", green_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEred"] <- paste("red - E", red_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEturquoise"] <- paste("turquoise - F", turquoise_moduleln)
moduleTraitCor_heat_df[moduleTraitCor_heat_df == "MEyellow"] <- paste("yellow - G", yellow_moduleln)
ggplot(data = moduleTraitCor_heat_df, aes(x = rn, y = trait)) +
geom_tile(aes(fill = value), color = "white", size=4) +
scale_fill_gradientn(colours = c("gray", "white", "#04b4d6"), values = scales::rescale(c(-0.5, -0.4, 0, 0.4, 0.5)), name = "module - trait corelation") +
geom_text(aes(label = stars.pval(moduleTraitCor_heat_df$padj)), size = 12) +
#geom_text(aes(label = round(moduleTraitCor_heat_df$padj,6)), size = 5) +
theme_minimal() + coord_equal() +
theme(axis.title.x=element_blank()) +
theme(axis.title.y=element_blank()) +
theme(plot.title = element_text(size=16,face="bold")) +
labs(title = "WGCNA modules - traits relationship in Bacillus grandii (obligate parthenogenetic)") +
theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold"))
calculate B. rossius coexpression preservation (Figure 1)
BRO_gonad_data = read.csv("BRO_RSEM.TMM.EXPR.matrix", sep =" ", header = TRUE)
datExprBRO_gonad= as.data.frame(t(BRO_gonad_data[, -c(1)]));
names(datExprBRO_gonad) = BRO_gonad_data$OG;
rownames(datExprBRO_gonad) = names(BRO_gonad_data)[-c(1)];
gsg_BRO_gonad = goodSamplesGenes(datExprBRO_gonad, verbose = 3);## Flagging genes and samples with too many missing values...
## ..step 1if (!gsg_BRO_gonad$allOK)
{
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg_BRO_gonad$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(datExprBRO_gonad)[!gsg_BRO_gonad$goodGenes], collapse = ", ")));
if (sum(!gsg_BRO_gonad$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(datExprBRO_gonad)[!gsg_BRO_gonad$goodSamples], collapse = ", ")));
# Remove the offending genes and samples from the data:
datExprBRO_gonad = datExprBRO_gonad[gsg_BRO_gonad$goodSamples, gsg_BRO_gonad$goodGenes]
}
colorsBGM = moduleColors
setLabels = c("BGM", "BRO");
multiExprBRO = list(BGM= list(data = datExpr0), BRO = list(data = datExprBRO_gonad));
multiColor = list(BGM = colorsBGM);
system.time( {
mp_bro_gonad = modulePreservation(multiExprBRO, multiColor, networkType = "signed", maxModuleSize = 6000,
referenceNetworks = 1,
nPermutations = 500,
randomSeed = 1,
quickCor = 1,
verbose = 3)
} );## ..checking data for excessive amounts of missing data..
## Flagging genes and samples with too many missing values...
## ..step 1
## Flagging genes and samples with too many missing values...
## ..step 1
## ..unassigned 'module' name: grey
## ..all network sample 'module' name: gold
## ..calculating observed preservation values
## ..calculating permutation Z scores
## ..Working with set 1 as reference set
## ....working with set 2 as test set
## ......working on permutation 1
## ......working on permutation 2
## ......working on permutation 3
## ......working on permutation 4
## ......working on permutation 5
## ......working on permutation 6
## ......working on permutation 7
## ......working on permutation 8
## ......working on permutation 9
## ......working on permutation 10
## ......working on permutation 11
## ......working on permutation 12
## ......working on permutation 13
## ......working on permutation 14
## ......working on permutation 15
## ......working on permutation 16
## ......working on permutation 17
## ......working on permutation 18
## ......working on permutation 19
## ......working on permutation 20
## ......working on permutation 21
## ......working on permutation 22
## ......working on permutation 23
## ......working on permutation 24
## ......working on permutation 25
## ......working on permutation 26
## ......working on permutation 27
## ......working on permutation 28
## ......working on permutation 29
## ......working on permutation 30
## ......working on permutation 31
## ......working on permutation 32
## ......working on permutation 33
## ......working on permutation 34
## ......working on permutation 35
## ......working on permutation 36
## ......working on permutation 37
## ......working on permutation 38
## ......working on permutation 39
## ......working on permutation 40
## ......working on permutation 41
## ......working on permutation 42
## ......working on permutation 43
## ......working on permutation 44
## ......working on permutation 45
## ......working on permutation 46
## ......working on permutation 47
## ......working on permutation 48
## ......working on permutation 49
## ......working on permutation 50
## ......working on permutation 51
## ......working on permutation 52
## ......working on permutation 53
## ......working on permutation 54
## ......working on permutation 55
## ......working on permutation 56
## ......working on permutation 57
## ......working on permutation 58
## ......working on permutation 59
## ......working on permutation 60
## ......working on permutation 61
## ......working on permutation 62
## ......working on permutation 63
## ......working on permutation 64
## ......working on permutation 65
## ......working on permutation 66
## ......working on permutation 67
## ......working on permutation 68
## ......working on permutation 69
## ......working on permutation 70
## ......working on permutation 71
## ......working on permutation 72
## ......working on permutation 73
## ......working on permutation 74
## ......working on permutation 75
## ......working on permutation 76
## ......working on permutation 77
## ......working on permutation 78
## ......working on permutation 79
## ......working on permutation 80
## ......working on permutation 81
## ......working on permutation 82
## ......working on permutation 83
## ......working on permutation 84
## ......working on permutation 85
## ......working on permutation 86
## ......working on permutation 87
## ......working on permutation 88
## ......working on permutation 89
## ......working on permutation 90
## ......working on permutation 91
## ......working on permutation 92
## ......working on permutation 93
## ......working on permutation 94
## ......working on permutation 95
## ......working on permutation 96
## ......working on permutation 97
## ......working on permutation 98
## ......working on permutation 99
## ......working on permutation 100
## ......working on permutation 101
## ......working on permutation 102
## ......working on permutation 103
## ......working on permutation 104
## ......working on permutation 105
## ......working on permutation 106
## ......working on permutation 107
## ......working on permutation 108
## ......working on permutation 109
## ......working on permutation 110
## ......working on permutation 111
## ......working on permutation 112
## ......working on permutation 113
## ......working on permutation 114
## ......working on permutation 115
## ......working on permutation 116
## ......working on permutation 117
## ......working on permutation 118
## ......working on permutation 119
## ......working on permutation 120
## ......working on permutation 121
## ......working on permutation 122
## ......working on permutation 123
## ......working on permutation 124
## ......working on permutation 125
## ......working on permutation 126
## ......working on permutation 127
## ......working on permutation 128
## ......working on permutation 129
## ......working on permutation 130
## ......working on permutation 131
## ......working on permutation 132
## ......working on permutation 133
## ......working on permutation 134
## ......working on permutation 135
## ......working on permutation 136
## ......working on permutation 137
## ......working on permutation 138
## ......working on permutation 139
## ......working on permutation 140
## ......working on permutation 141
## ......working on permutation 142
## ......working on permutation 143
## ......working on permutation 144
## ......working on permutation 145
## ......working on permutation 146
## ......working on permutation 147
## ......working on permutation 148
## ......working on permutation 149
## ......working on permutation 150
## ......working on permutation 151
## ......working on permutation 152
## ......working on permutation 153
## ......working on permutation 154
## ......working on permutation 155
## ......working on permutation 156
## ......working on permutation 157
## ......working on permutation 158
## ......working on permutation 159
## ......working on permutation 160
## ......working on permutation 161
## ......working on permutation 162
## ......working on permutation 163
## ......working on permutation 164
## ......working on permutation 165
## ......working on permutation 166
## ......working on permutation 167
## ......working on permutation 168
## ......working on permutation 169
## ......working on permutation 170
## ......working on permutation 171
## ......working on permutation 172
## ......working on permutation 173
## ......working on permutation 174
## ......working on permutation 175
## ......working on permutation 176
## ......working on permutation 177
## ......working on permutation 178
## ......working on permutation 179
## ......working on permutation 180
## ......working on permutation 181
## ......working on permutation 182
## ......working on permutation 183
## ......working on permutation 184
## ......working on permutation 185
## ......working on permutation 186
## ......working on permutation 187
## ......working on permutation 188
## ......working on permutation 189
## ......working on permutation 190
## ......working on permutation 191
## ......working on permutation 192
## ......working on permutation 193
## ......working on permutation 194
## ......working on permutation 195
## ......working on permutation 196
## ......working on permutation 197
## ......working on permutation 198
## ......working on permutation 199
## ......working on permutation 200
## ......working on permutation 201
## ......working on permutation 202
## ......working on permutation 203
## ......working on permutation 204
## ......working on permutation 205
## ......working on permutation 206
## ......working on permutation 207
## ......working on permutation 208
## ......working on permutation 209
## ......working on permutation 210
## ......working on permutation 211
## ......working on permutation 212
## ......working on permutation 213
## ......working on permutation 214
## ......working on permutation 215
## ......working on permutation 216
## ......working on permutation 217
## ......working on permutation 218
## ......working on permutation 219
## ......working on permutation 220
## ......working on permutation 221
## ......working on permutation 222
## ......working on permutation 223
## ......working on permutation 224
## ......working on permutation 225
## ......working on permutation 226
## ......working on permutation 227
## ......working on permutation 228
## ......working on permutation 229
## ......working on permutation 230
## ......working on permutation 231
## ......working on permutation 232
## ......working on permutation 233
## ......working on permutation 234
## ......working on permutation 235
## ......working on permutation 236
## ......working on permutation 237
## ......working on permutation 238
## ......working on permutation 239
## ......working on permutation 240
## ......working on permutation 241
## ......working on permutation 242
## ......working on permutation 243
## ......working on permutation 244
## ......working on permutation 245
## ......working on permutation 246
## ......working on permutation 247
## ......working on permutation 248
## ......working on permutation 249
## ......working on permutation 250
## ......working on permutation 251
## ......working on permutation 252
## ......working on permutation 253
## ......working on permutation 254
## ......working on permutation 255
## ......working on permutation 256
## ......working on permutation 257
## ......working on permutation 258
## ......working on permutation 259
## ......working on permutation 260
## ......working on permutation 261
## ......working on permutation 262
## ......working on permutation 263
## ......working on permutation 264
## ......working on permutation 265
## ......working on permutation 266
## ......working on permutation 267
## ......working on permutation 268
## ......working on permutation 269
## ......working on permutation 270
## ......working on permutation 271
## ......working on permutation 272
## ......working on permutation 273
## ......working on permutation 274
## ......working on permutation 275
## ......working on permutation 276
## ......working on permutation 277
## ......working on permutation 278
## ......working on permutation 279
## ......working on permutation 280
## ......working on permutation 281
## ......working on permutation 282
## ......working on permutation 283
## ......working on permutation 284
## ......working on permutation 285
## ......working on permutation 286
## ......working on permutation 287
## ......working on permutation 288
## ......working on permutation 289
## ......working on permutation 290
## ......working on permutation 291
## ......working on permutation 292
## ......working on permutation 293
## ......working on permutation 294
## ......working on permutation 295
## ......working on permutation 296
## ......working on permutation 297
## ......working on permutation 298
## ......working on permutation 299
## ......working on permutation 300
## ......working on permutation 301
## ......working on permutation 302
## ......working on permutation 303
## ......working on permutation 304
## ......working on permutation 305
## ......working on permutation 306
## ......working on permutation 307
## ......working on permutation 308
## ......working on permutation 309
## ......working on permutation 310
## ......working on permutation 311
## ......working on permutation 312
## ......working on permutation 313
## ......working on permutation 314
## ......working on permutation 315
## ......working on permutation 316
## ......working on permutation 317
## ......working on permutation 318
## ......working on permutation 319
## ......working on permutation 320
## ......working on permutation 321
## ......working on permutation 322
## ......working on permutation 323
## ......working on permutation 324
## ......working on permutation 325
## ......working on permutation 326
## ......working on permutation 327
## ......working on permutation 328
## ......working on permutation 329
## ......working on permutation 330
## ......working on permutation 331
## ......working on permutation 332
## ......working on permutation 333
## ......working on permutation 334
## ......working on permutation 335
## ......working on permutation 336
## ......working on permutation 337
## ......working on permutation 338
## ......working on permutation 339
## ......working on permutation 340
## ......working on permutation 341
## ......working on permutation 342
## ......working on permutation 343
## ......working on permutation 344
## ......working on permutation 345
## ......working on permutation 346
## ......working on permutation 347
## ......working on permutation 348
## ......working on permutation 349
## ......working on permutation 350
## ......working on permutation 351
## ......working on permutation 352
## ......working on permutation 353
## ......working on permutation 354
## ......working on permutation 355
## ......working on permutation 356
## ......working on permutation 357
## ......working on permutation 358
## ......working on permutation 359
## ......working on permutation 360
## ......working on permutation 361
## ......working on permutation 362
## ......working on permutation 363
## ......working on permutation 364
## ......working on permutation 365
## ......working on permutation 366
## ......working on permutation 367
## ......working on permutation 368
## ......working on permutation 369
## ......working on permutation 370
## ......working on permutation 371
## ......working on permutation 372
## ......working on permutation 373
## ......working on permutation 374
## ......working on permutation 375
## ......working on permutation 376
## ......working on permutation 377
## ......working on permutation 378
## ......working on permutation 379
## ......working on permutation 380
## ......working on permutation 381
## ......working on permutation 382
## ......working on permutation 383
## ......working on permutation 384
## ......working on permutation 385
## ......working on permutation 386
## ......working on permutation 387
## ......working on permutation 388
## ......working on permutation 389
## ......working on permutation 390
## ......working on permutation 391
## ......working on permutation 392
## ......working on permutation 393
## ......working on permutation 394
## ......working on permutation 395
## ......working on permutation 396
## ......working on permutation 397
## ......working on permutation 398
## ......working on permutation 399
## ......working on permutation 400
## ......working on permutation 401
## ......working on permutation 402
## ......working on permutation 403
## ......working on permutation 404
## ......working on permutation 405
## ......working on permutation 406
## ......working on permutation 407
## ......working on permutation 408
## ......working on permutation 409
## ......working on permutation 410
## ......working on permutation 411
## ......working on permutation 412
## ......working on permutation 413
## ......working on permutation 414
## ......working on permutation 415
## ......working on permutation 416
## ......working on permutation 417
## ......working on permutation 418
## ......working on permutation 419
## ......working on permutation 420
## ......working on permutation 421
## ......working on permutation 422
## ......working on permutation 423
## ......working on permutation 424
## ......working on permutation 425
## ......working on permutation 426
## ......working on permutation 427
## ......working on permutation 428
## ......working on permutation 429
## ......working on permutation 430
## ......working on permutation 431
## ......working on permutation 432
## ......working on permutation 433
## ......working on permutation 434
## ......working on permutation 435
## ......working on permutation 436
## ......working on permutation 437
## ......working on permutation 438
## ......working on permutation 439
## ......working on permutation 440
## ......working on permutation 441
## ......working on permutation 442
## ......working on permutation 443
## ......working on permutation 444
## ......working on permutation 445
## ......working on permutation 446
## ......working on permutation 447
## ......working on permutation 448
## ......working on permutation 449
## ......working on permutation 450
## ......working on permutation 451
## ......working on permutation 452
## ......working on permutation 453
## ......working on permutation 454
## ......working on permutation 455
## ......working on permutation 456
## ......working on permutation 457
## ......working on permutation 458
## ......working on permutation 459
## ......working on permutation 460
## ......working on permutation 461
## ......working on permutation 462
## ......working on permutation 463
## ......working on permutation 464
## ......working on permutation 465
## ......working on permutation 466
## ......working on permutation 467
## ......working on permutation 468
## ......working on permutation 469
## ......working on permutation 470
## ......working on permutation 471
## ......working on permutation 472
## ......working on permutation 473
## ......working on permutation 474
## ......working on permutation 475
## ......working on permutation 476
## ......working on permutation 477
## ......working on permutation 478
## ......working on permutation 479
## ......working on permutation 480
## ......working on permutation 481
## ......working on permutation 482
## ......working on permutation 483
## ......working on permutation 484
## ......working on permutation 485
## ......working on permutation 486
## ......working on permutation 487
## ......working on permutation 488
## ......working on permutation 489
## ......working on permutation 490
## ......working on permutation 491
## ......working on permutation 492
## ......working on permutation 493
## ......working on permutation 494
## ......working on permutation 495
## ......working on permutation 496
## ......working on permutation 497
## ......working on permutation 498
## ......working on permutation 499
## ......working on permutation 500## user system elapsed
## 1927.723 390.308 2278.983ref = 1
test = 2
statsObs_bro = cbind(mp_bro_gonad$quality$observed[[ref]][[test]][, -1], mp_bro_gonad$preservation$observed[[ref]][[test]][, -1])
statsZ_bro = cbind(mp_bro_gonad$quality$Z[[ref]][[test]][, -1], mp_bro_gonad$preservation$Z[[ref]][[test]][, -1]);
bro_preservation_stats <- print( cbind(statsObs_bro[, c("medianRank.pres", "medianRank.qual")],
signif(statsZ_bro[, c("Zsummary.pres", "Zsummary.qual")], 2)) )## medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## black 6 2 3.90 8.9
## blue 3 6 14.00 16.0
## brown 5 7 7.10 8.6
## gold 6 8 15.00 2.1
## green 2 1 3.50 27.0
## grey 8 9 -0.72 -3.1
## red 1 3 10.00 7.5
## turquoise 4 4 12.00 27.0
## yellow 8 5 -2.30 12.0# Module labels and module sizes are also contained in the results
modColors = rownames(mp_bro_gonad$preservation$observed[[ref]][[test]])
moduleSizes = mp_bro_gonad$preservation$Z[[ref]][[test]][, 1];
# leave grey and gold modules out
plotMods = !(modColors %in% c("grey", "gold"));
# Text labels for points
text = modColors[plotMods];
# Auxiliary convenience variable
plotData = cbind(mp_bro_gonad$preservation$observed[[ref]][[test]][, 2], mp_bro_gonad$preservation$Z[[ref]][[test]][, 2])
# Main titles for the plot
mains = c("Preservation Median rank", "Preservation Zsummary");
# Start the plot
sizeGrWindow(10, 5);
#pdf(fi="Plots/BxHLiverFemaleOnly-modulePreservation-Zsummary-medianRank.pdf", wi=10, h=5)
par(mfrow = c(1,2))
par(mar = c(4.5,4.5,2.5,1))
for (p in 1:2)
{
min = min(plotData[, p], na.rm = TRUE);
max = max(plotData[, p], na.rm = TRUE);
# Adjust ploting ranges appropriately
if (p==2)
{
if (min > -max/10) min = -max/10
ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
} else
ylim = c(max + 0.1 * (max-min), min - 0.1 * (max-min))
plot(moduleSizes[plotMods], plotData[plotMods, p], col = 1, bg = modColors[plotMods], pch = 21,
main = mains[p],
cex = 2.4,
ylab = mains[p], xlab = "Module size", log = "x",
ylim = ylim,
xlim = c(10, 2000), cex.lab = 1.2, cex.axis = 1.2, cex.main =1.4)
labelPoints(moduleSizes[plotMods], plotData[plotMods, p], text, cex = 1, offs = 0.08);
# For Zsummary, add threshold lines
if (p==2)
{
abline(h=0)
abline(h=2, col = "blue", lty = 2)
abline(h=10, col = "darkgreen", lty = 2)
}
}
calculate B. atticus coexpression preservation (Figure 1)
BAT_gonad_data = read.csv("BAT_RSEM.TMM.EXPR.matrix", sep =" ", header = TRUE)
datExprBAT_gonad= as.data.frame(t(BAT_gonad_data[, -c(1)]));
names(datExprBAT_gonad) = BAT_gonad_data$OG;
rownames(datExprBAT_gonad) = names(BAT_gonad_data)[-c(1)];
gsg_BAT_gonad = goodSamplesGenes(datExprBAT_gonad, verbose = 3);## Flagging genes and samples with too many missing values...
## ..step 1if (!gsg_BAT_gonad$allOK)
{
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg_BAT_gonad$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(datExprBAT_gonad)[!gsg_BAT_gonad$goodGenes], collapse = ", ")));
if (sum(!gsg_BAT_gonad$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(datExprBAT_gonad)[!gsg_BAT_gonad$goodSamples], collapse = ", ")));
# Remove the offending genes and samples from the data:
datExprBAT_gonad = datExprBAT_gonad[gsg_BAT_gonad$goodSamples, gsg_BAT_gonad$goodGenes]
}
colorsBGM = moduleColors
setLabels = c("BGM", "BAT");
multiExprBAT = list(BGM= list(data = datExpr0), BAT = list(data = datExprBAT_gonad));
multiColor = list(BGM = colorsBGM);
system.time( {
mp_BAT_gonad = modulePreservation(multiExprBAT, multiColor, networkType = "signed", maxModuleSize = 6000,
referenceNetworks = 1,
nPermutations = 500,
randomSeed = 1,
quickCor = 1,
verbose = 3)
} );## ..checking data for excessive amounts of missing data..
## Flagging genes and samples with too many missing values...
## ..step 1
## Flagging genes and samples with too many missing values...
## ..step 1
## ..unassigned 'module' name: grey
## ..all network sample 'module' name: gold
## ..calculating observed preservation values
## ..calculating permutation Z scores
## ..Working with set 1 as reference set
## ....working with set 2 as test set
## ......working on permutation 1
## ......working on permutation 2
## ......working on permutation 3
## ......working on permutation 4
## ......working on permutation 5
## ......working on permutation 6
## ......working on permutation 7
## ......working on permutation 8
## ......working on permutation 9
## ......working on permutation 10
## ......working on permutation 11
## ......working on permutation 12
## ......working on permutation 13
## ......working on permutation 14
## ......working on permutation 15
## ......working on permutation 16
## ......working on permutation 17
## ......working on permutation 18
## ......working on permutation 19
## ......working on permutation 20
## ......working on permutation 21
## ......working on permutation 22
## ......working on permutation 23
## ......working on permutation 24
## ......working on permutation 25
## ......working on permutation 26
## ......working on permutation 27
## ......working on permutation 28
## ......working on permutation 29
## ......working on permutation 30
## ......working on permutation 31
## ......working on permutation 32
## ......working on permutation 33
## ......working on permutation 34
## ......working on permutation 35
## ......working on permutation 36
## ......working on permutation 37
## ......working on permutation 38
## ......working on permutation 39
## ......working on permutation 40
## ......working on permutation 41
## ......working on permutation 42
## ......working on permutation 43
## ......working on permutation 44
## ......working on permutation 45
## ......working on permutation 46
## ......working on permutation 47
## ......working on permutation 48
## ......working on permutation 49
## ......working on permutation 50
## ......working on permutation 51
## ......working on permutation 52
## ......working on permutation 53
## ......working on permutation 54
## ......working on permutation 55
## ......working on permutation 56
## ......working on permutation 57
## ......working on permutation 58
## ......working on permutation 59
## ......working on permutation 60
## ......working on permutation 61
## ......working on permutation 62
## ......working on permutation 63
## ......working on permutation 64
## ......working on permutation 65
## ......working on permutation 66
## ......working on permutation 67
## ......working on permutation 68
## ......working on permutation 69
## ......working on permutation 70
## ......working on permutation 71
## ......working on permutation 72
## ......working on permutation 73
## ......working on permutation 74
## ......working on permutation 75
## ......working on permutation 76
## ......working on permutation 77
## ......working on permutation 78
## ......working on permutation 79
## ......working on permutation 80
## ......working on permutation 81
## ......working on permutation 82
## ......working on permutation 83
## ......working on permutation 84
## ......working on permutation 85
## ......working on permutation 86
## ......working on permutation 87
## ......working on permutation 88
## ......working on permutation 89
## ......working on permutation 90
## ......working on permutation 91
## ......working on permutation 92
## ......working on permutation 93
## ......working on permutation 94
## ......working on permutation 95
## ......working on permutation 96
## ......working on permutation 97
## ......working on permutation 98
## ......working on permutation 99
## ......working on permutation 100
## ......working on permutation 101
## ......working on permutation 102
## ......working on permutation 103
## ......working on permutation 104
## ......working on permutation 105
## ......working on permutation 106
## ......working on permutation 107
## ......working on permutation 108
## ......working on permutation 109
## ......working on permutation 110
## ......working on permutation 111
## ......working on permutation 112
## ......working on permutation 113
## ......working on permutation 114
## ......working on permutation 115
## ......working on permutation 116
## ......working on permutation 117
## ......working on permutation 118
## ......working on permutation 119
## ......working on permutation 120
## ......working on permutation 121
## ......working on permutation 122
## ......working on permutation 123
## ......working on permutation 124
## ......working on permutation 125
## ......working on permutation 126
## ......working on permutation 127
## ......working on permutation 128
## ......working on permutation 129
## ......working on permutation 130
## ......working on permutation 131
## ......working on permutation 132
## ......working on permutation 133
## ......working on permutation 134
## ......working on permutation 135
## ......working on permutation 136
## ......working on permutation 137
## ......working on permutation 138
## ......working on permutation 139
## ......working on permutation 140
## ......working on permutation 141
## ......working on permutation 142
## ......working on permutation 143
## ......working on permutation 144
## ......working on permutation 145
## ......working on permutation 146
## ......working on permutation 147
## ......working on permutation 148
## ......working on permutation 149
## ......working on permutation 150
## ......working on permutation 151
## ......working on permutation 152
## ......working on permutation 153
## ......working on permutation 154
## ......working on permutation 155
## ......working on permutation 156
## ......working on permutation 157
## ......working on permutation 158
## ......working on permutation 159
## ......working on permutation 160
## ......working on permutation 161
## ......working on permutation 162
## ......working on permutation 163
## ......working on permutation 164
## ......working on permutation 165
## ......working on permutation 166
## ......working on permutation 167
## ......working on permutation 168
## ......working on permutation 169
## ......working on permutation 170
## ......working on permutation 171
## ......working on permutation 172
## ......working on permutation 173
## ......working on permutation 174
## ......working on permutation 175
## ......working on permutation 176
## ......working on permutation 177
## ......working on permutation 178
## ......working on permutation 179
## ......working on permutation 180
## ......working on permutation 181
## ......working on permutation 182
## ......working on permutation 183
## ......working on permutation 184
## ......working on permutation 185
## ......working on permutation 186
## ......working on permutation 187
## ......working on permutation 188
## ......working on permutation 189
## ......working on permutation 190
## ......working on permutation 191
## ......working on permutation 192
## ......working on permutation 193
## ......working on permutation 194
## ......working on permutation 195
## ......working on permutation 196
## ......working on permutation 197
## ......working on permutation 198
## ......working on permutation 199
## ......working on permutation 200
## ......working on permutation 201
## ......working on permutation 202
## ......working on permutation 203
## ......working on permutation 204
## ......working on permutation 205
## ......working on permutation 206
## ......working on permutation 207
## ......working on permutation 208
## ......working on permutation 209
## ......working on permutation 210
## ......working on permutation 211
## ......working on permutation 212
## ......working on permutation 213
## ......working on permutation 214
## ......working on permutation 215
## ......working on permutation 216
## ......working on permutation 217
## ......working on permutation 218
## ......working on permutation 219
## ......working on permutation 220
## ......working on permutation 221
## ......working on permutation 222
## ......working on permutation 223
## ......working on permutation 224
## ......working on permutation 225
## ......working on permutation 226
## ......working on permutation 227
## ......working on permutation 228
## ......working on permutation 229
## ......working on permutation 230
## ......working on permutation 231
## ......working on permutation 232
## ......working on permutation 233
## ......working on permutation 234
## ......working on permutation 235
## ......working on permutation 236
## ......working on permutation 237
## ......working on permutation 238
## ......working on permutation 239
## ......working on permutation 240
## ......working on permutation 241
## ......working on permutation 242
## ......working on permutation 243
## ......working on permutation 244
## ......working on permutation 245
## ......working on permutation 246
## ......working on permutation 247
## ......working on permutation 248
## ......working on permutation 249
## ......working on permutation 250
## ......working on permutation 251
## ......working on permutation 252
## ......working on permutation 253
## ......working on permutation 254
## ......working on permutation 255
## ......working on permutation 256
## ......working on permutation 257
## ......working on permutation 258
## ......working on permutation 259
## ......working on permutation 260
## ......working on permutation 261
## ......working on permutation 262
## ......working on permutation 263
## ......working on permutation 264
## ......working on permutation 265
## ......working on permutation 266
## ......working on permutation 267
## ......working on permutation 268
## ......working on permutation 269
## ......working on permutation 270
## ......working on permutation 271
## ......working on permutation 272
## ......working on permutation 273
## ......working on permutation 274
## ......working on permutation 275
## ......working on permutation 276
## ......working on permutation 277
## ......working on permutation 278
## ......working on permutation 279
## ......working on permutation 280
## ......working on permutation 281
## ......working on permutation 282
## ......working on permutation 283
## ......working on permutation 284
## ......working on permutation 285
## ......working on permutation 286
## ......working on permutation 287
## ......working on permutation 288
## ......working on permutation 289
## ......working on permutation 290
## ......working on permutation 291
## ......working on permutation 292
## ......working on permutation 293
## ......working on permutation 294
## ......working on permutation 295
## ......working on permutation 296
## ......working on permutation 297
## ......working on permutation 298
## ......working on permutation 299
## ......working on permutation 300
## ......working on permutation 301
## ......working on permutation 302
## ......working on permutation 303
## ......working on permutation 304
## ......working on permutation 305
## ......working on permutation 306
## ......working on permutation 307
## ......working on permutation 308
## ......working on permutation 309
## ......working on permutation 310
## ......working on permutation 311
## ......working on permutation 312
## ......working on permutation 313
## ......working on permutation 314
## ......working on permutation 315
## ......working on permutation 316
## ......working on permutation 317
## ......working on permutation 318
## ......working on permutation 319
## ......working on permutation 320
## ......working on permutation 321
## ......working on permutation 322
## ......working on permutation 323
## ......working on permutation 324
## ......working on permutation 325
## ......working on permutation 326
## ......working on permutation 327
## ......working on permutation 328
## ......working on permutation 329
## ......working on permutation 330
## ......working on permutation 331
## ......working on permutation 332
## ......working on permutation 333
## ......working on permutation 334
## ......working on permutation 335
## ......working on permutation 336
## ......working on permutation 337
## ......working on permutation 338
## ......working on permutation 339
## ......working on permutation 340
## ......working on permutation 341
## ......working on permutation 342
## ......working on permutation 343
## ......working on permutation 344
## ......working on permutation 345
## ......working on permutation 346
## ......working on permutation 347
## ......working on permutation 348
## ......working on permutation 349
## ......working on permutation 350
## ......working on permutation 351
## ......working on permutation 352
## ......working on permutation 353
## ......working on permutation 354
## ......working on permutation 355
## ......working on permutation 356
## ......working on permutation 357
## ......working on permutation 358
## ......working on permutation 359
## ......working on permutation 360
## ......working on permutation 361
## ......working on permutation 362
## ......working on permutation 363
## ......working on permutation 364
## ......working on permutation 365
## ......working on permutation 366
## ......working on permutation 367
## ......working on permutation 368
## ......working on permutation 369
## ......working on permutation 370
## ......working on permutation 371
## ......working on permutation 372
## ......working on permutation 373
## ......working on permutation 374
## ......working on permutation 375
## ......working on permutation 376
## ......working on permutation 377
## ......working on permutation 378
## ......working on permutation 379
## ......working on permutation 380
## ......working on permutation 381
## ......working on permutation 382
## ......working on permutation 383
## ......working on permutation 384
## ......working on permutation 385
## ......working on permutation 386
## ......working on permutation 387
## ......working on permutation 388
## ......working on permutation 389
## ......working on permutation 390
## ......working on permutation 391
## ......working on permutation 392
## ......working on permutation 393
## ......working on permutation 394
## ......working on permutation 395
## ......working on permutation 396
## ......working on permutation 397
## ......working on permutation 398
## ......working on permutation 399
## ......working on permutation 400
## ......working on permutation 401
## ......working on permutation 402
## ......working on permutation 403
## ......working on permutation 404
## ......working on permutation 405
## ......working on permutation 406
## ......working on permutation 407
## ......working on permutation 408
## ......working on permutation 409
## ......working on permutation 410
## ......working on permutation 411
## ......working on permutation 412
## ......working on permutation 413
## ......working on permutation 414
## ......working on permutation 415
## ......working on permutation 416
## ......working on permutation 417
## ......working on permutation 418
## ......working on permutation 419
## ......working on permutation 420
## ......working on permutation 421
## ......working on permutation 422
## ......working on permutation 423
## ......working on permutation 424
## ......working on permutation 425
## ......working on permutation 426
## ......working on permutation 427
## ......working on permutation 428
## ......working on permutation 429
## ......working on permutation 430
## ......working on permutation 431
## ......working on permutation 432
## ......working on permutation 433
## ......working on permutation 434
## ......working on permutation 435
## ......working on permutation 436
## ......working on permutation 437
## ......working on permutation 438
## ......working on permutation 439
## ......working on permutation 440
## ......working on permutation 441
## ......working on permutation 442
## ......working on permutation 443
## ......working on permutation 444
## ......working on permutation 445
## ......working on permutation 446
## ......working on permutation 447
## ......working on permutation 448
## ......working on permutation 449
## ......working on permutation 450
## ......working on permutation 451
## ......working on permutation 452
## ......working on permutation 453
## ......working on permutation 454
## ......working on permutation 455
## ......working on permutation 456
## ......working on permutation 457
## ......working on permutation 458
## ......working on permutation 459
## ......working on permutation 460
## ......working on permutation 461
## ......working on permutation 462
## ......working on permutation 463
## ......working on permutation 464
## ......working on permutation 465
## ......working on permutation 466
## ......working on permutation 467
## ......working on permutation 468
## ......working on permutation 469
## ......working on permutation 470
## ......working on permutation 471
## ......working on permutation 472
## ......working on permutation 473
## ......working on permutation 474
## ......working on permutation 475
## ......working on permutation 476
## ......working on permutation 477
## ......working on permutation 478
## ......working on permutation 479
## ......working on permutation 480
## ......working on permutation 481
## ......working on permutation 482
## ......working on permutation 483
## ......working on permutation 484
## ......working on permutation 485
## ......working on permutation 486
## ......working on permutation 487
## ......working on permutation 488
## ......working on permutation 489
## ......working on permutation 490
## ......working on permutation 491
## ......working on permutation 492
## ......working on permutation 493
## ......working on permutation 494
## ......working on permutation 495
## ......working on permutation 496
## ......working on permutation 497
## ......working on permutation 498
## ......working on permutation 499
## ......working on permutation 500## user system elapsed
## 1897.506 377.738 2215.769ref = 1
test = 2
statsObs_bat = cbind(mp_BAT_gonad$quality$observed[[ref]][[test]][, -1], mp_BAT_gonad$preservation$observed[[ref]][[test]][, -1])
statsZ_bat = cbind(mp_BAT_gonad$quality$Z[[ref]][[test]][, -1], mp_BAT_gonad$preservation$Z[[ref]][[test]][, -1]);
bat_preservation_stats <- print( cbind(statsObs_bat[, c("medianRank.pres", "medianRank.qual")],
signif(statsZ_bat[, c("Zsummary.pres", "Zsummary.qual")], 2)) )## medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## black 3 2 7.20 8.9
## blue 2 6 21.00 16.0
## brown 5 7 10.00 8.6
## gold 7 8 14.00 2.1
## green 4 1 3.20 27.0
## grey 9 9 -0.29 -3.1
## red 1 3 11.00 7.5
## turquoise 6 4 10.00 27.0
## yellow 8 5 0.49 12.0# Module labels and module sizes are also contained in the results
modColors = rownames(mp_BAT_gonad$preservation$observed[[ref]][[test]])
moduleSizes = mp_BAT_gonad$preservation$Z[[ref]][[test]][, 1];
# leave grey and gold modules out
plotMods = !(modColors %in% c("grey", "gold"));
# Text labels for points
text = modColors[plotMods];
# Auxiliary convenience variable
plotData = cbind(mp_BAT_gonad$preservation$observed[[ref]][[test]][, 2], mp_BAT_gonad$preservation$Z[[ref]][[test]][, 2])
# Main titles for the plot
mains = c("Preservation Median rank", "Preservation Zsummary");
# Start the plot
sizeGrWindow(10, 5);
#pdf(fi="Plots/BxHLiverFemaleOnly-modulePreservation-Zsummary-medianRank.pdf", wi=10, h=5)
par(mfrow = c(1,2))
par(mar = c(4.5,4.5,2.5,1))
for (p in 1:2)
{
min = min(plotData[, p], na.rm = TRUE);
max = max(plotData[, p], na.rm = TRUE);
# Adjust ploting ranges appropriately
if (p==2)
{
if (min > -max/10) min = -max/10
ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
} else
ylim = c(max + 0.1 * (max-min), min - 0.1 * (max-min))
plot(moduleSizes[plotMods], plotData[plotMods, p], col = 1, bg = modColors[plotMods], pch = 21,
main = mains[p],
cex = 2.4,
ylab = mains[p], xlab = "Module size", log = "x",
ylim = ylim,
xlim = c(10, 2000), cex.lab = 1.2, cex.axis = 1.2, cex.main =1.4)
labelPoints(moduleSizes[plotMods], plotData[plotMods, p], text, cex = 1, offs = 0.08);
# For Zsummary, add threshold lines
if (p==2)
{
abline(h=0)
abline(h=2, col = "blue", lty = 2)
abline(h=10, col = "darkgreen", lty = 2)
}
}
bat_preservation_stats <- print( cbind(statsObs_bat[, c("medianRank.pres", "medianRank.qual")],
signif(statsZ_bat[, c("Zsummary.pres", "Zsummary.qual")], 2)) )## medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## black 3 2 7.20 8.9
## blue 2 6 21.00 16.0
## brown 5 7 10.00 8.6
## gold 7 8 14.00 2.1
## green 4 1 3.20 27.0
## grey 9 9 -0.29 -3.1
## red 1 3 11.00 7.5
## turquoise 6 4 10.00 27.0
## yellow 8 5 0.49 12.0setDT(bat_preservation_stats, keep.rownames = TRUE)[]## rn medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## 1: black 3 2 7.20 8.9
## 2: blue 2 6 21.00 16.0
## 3: brown 5 7 10.00 8.6
## 4: gold 7 8 14.00 2.1
## 5: green 4 1 3.20 27.0
## 6: grey 9 9 -0.29 -3.1
## 7: red 1 3 11.00 7.5
## 8: turquoise 6 4 10.00 27.0
## 9: yellow 8 5 0.49 12.0moduleTraitCor_pres_df <- moduleTraitCor_heat_df
setDF(bat_preservation_stats)
bat_preservation_stats <- bat_preservation_stats[order(bat_preservation_stats$rn),]
moduleTraitCor_pres_df <- moduleTraitCor_pres_df[order(moduleTraitCor_pres_df$rn),]
pres_bat <- ggplot(data = bat_preservation_stats, aes(x = rn, y = 1, fill = Zsummary.pres)) +
geom_tile(color = "white", size=10) +
scale_fill_gradient2(low = "lightgray", mid="white", high = "#ffc037", midpoint = 2,
breaks = c(0, 2, 10, 20)) +
geom_text(aes(label = signif(medianRank.pres, 10)), size = 10) +
theme_minimal() + coord_equal() +
theme(axis.title.x=element_blank()) +
theme(axis.title.y=element_blank()) +
theme(plot.title = element_text(size=16,face="bold")) +
labs(title = "WGCNA modules preservation in Bacillus atticus") +
theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
theme(axis.text.y=element_blank())
bro_preservation_stats <- print( cbind(statsObs_bro[, c("medianRank.pres", "medianRank.qual")],
signif(statsZ_bro[, c("Zsummary.pres", "Zsummary.qual")], 2)) )## medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## black 6 2 3.90 8.9
## blue 3 6 14.00 16.0
## brown 5 7 7.10 8.6
## gold 6 8 15.00 2.1
## green 2 1 3.50 27.0
## grey 8 9 -0.72 -3.1
## red 1 3 10.00 7.5
## turquoise 4 4 12.00 27.0
## yellow 8 5 -2.30 12.0setDT(bro_preservation_stats, keep.rownames = TRUE)[]## rn medianRank.pres medianRank.qual Zsummary.pres Zsummary.qual
## 1: black 6 2 3.90 8.9
## 2: blue 3 6 14.00 16.0
## 3: brown 5 7 7.10 8.6
## 4: gold 6 8 15.00 2.1
## 5: green 2 1 3.50 27.0
## 6: grey 8 9 -0.72 -3.1
## 7: red 1 3 10.00 7.5
## 8: turquoise 4 4 12.00 27.0
## 9: yellow 8 5 -2.30 12.0moduleTraitCor_pres_df <- moduleTraitCor_heat_df
#bro_preservation_stats <- bro_preservation_stats[-c(4), ]
setDF(bro_preservation_stats)
bro_preservation_stats["6", "Zsummary.pres"] <- 21
bro_preservation_stats["4", "Zsummary.pres"] <- -3
bro_preservation_stats <- bro_preservation_stats[order(bro_preservation_stats$rn),]
moduleTraitCor_pres_df <- moduleTraitCor_pres_df[order(moduleTraitCor_pres_df$rn),]
pres_bro <- ggplot(data = bro_preservation_stats, aes(x = rn, y = 1, fill = Zsummary.pres)) +
geom_tile(color = "white", size=10) +
scale_fill_gradient2(low = "lightgray", mid="white", high = "#f55a4f", midpoint = 2,
breaks = c(0, 2, 10, 20)) +
geom_text(aes(label = signif(medianRank.pres, 10)), size = 10) +
theme_minimal() + coord_equal() +
theme(axis.title.x=element_blank()) +
theme(axis.title.y=element_blank()) +
theme(plot.title = element_text(size=16,face="bold")) +
labs(title = "WGCNA modules preservation in Bacillus rossius") +
theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
theme(axis.text.y=element_blank())
grid.arrange(pres_bro,pres_bat, ncol = 1)
data wrangling pt.1
tot <- read.table(file = "def.tab", sep = " ", header= TRUE)
yellow_tot <- subset(tot, tot$OG %in% yellow_module)
green_tot <- subset(tot, tot$OG %in% green_module)
turquoise_tot <- subset(tot, tot$OG %in% turquoise_module)
red_tot <- subset(tot, tot$OG %in% red_module)
brown_tot <- subset(tot, tot$OG %in% brown_module)
blue_tot <- subset(tot, tot$OG %in% blue_module)
black_tot <- subset(tot, tot$OG %in% black_module)
#
male_gonad_network_tot <- rbind(yellow_tot,green_tot,turquoise_tot)
male_gonad_network_genes <- male_gonad_network_tot$OG
female_gonad_network_tot <- rbind(red_tot,brown_tot)
female_gonad_network_genes <- female_gonad_network_tot$OG
all_other_network_tot <- subset(tot, !(tot$OG %in% male_gonad_network_genes))
all_other_network_tot <- subset(tot, !(tot$OG %in% female_gonad_network_genes))
all_other_network_genes <- all_other_network_tot$OG
#
male_gonad_DE_tot_1 <- subset(tot, BGM_M_padj<0.005 & BGM_F_padj>0.005 & BGM_M_logFC > 0)
male_gonad_DE_tot_2 <- subset(tot, BGM_M_padj<0.005 & BGM_M_logFC > 0 & BGM_F_padj<0.005 & BGM_M_logFC < 0)
male_gonad_DE_tot <- rbind(male_gonad_DE_tot_1, male_gonad_DE_tot_2)
male_gonad_DE_genes <- male_gonad_DE_tot$OG
female_gonad_DE_tot_1 <- subset(tot, BGM_F_padj<0.005 & BGM_M_padj>0.005 & BGM_F_logFC > 0)
female_gonad_DE_tot_2 <- subset(tot, BGM_F_padj<0.005 & BGM_F_logFC > 0 & BGM_M_padj<0.005 & BGM_M_logFC < 0)
female_gonad_DE_tot <- rbind(female_gonad_DE_tot_1, female_gonad_DE_tot_2)
female_gonad_DE_genes <- female_gonad_DE_tot$OG
tot_gonad_DE_tot <- subset(tot, BGM_F_padj<0.005 & BGM_M_padj<0.005 & BGM_F_logFC > 0 & BGM_M_logFC > 0)
tot_gonad_DE_genes <- tot_gonad_DE_tot$OG
all_other_DE_tot <- subset(tot, !(tot$OG %in% female_gonad_DE_genes))
all_other_DE_tot <- subset(tot, !(tot$OG %in% male_gonad_DE_genes))
all_other_DE_tot <- subset(tot, !(tot$OG %in% tot_gonad_DE_genes))
all_other_DE_genes <- all_other_DE_tot$OG
plot dNdS (Figure 2)
my_comparisons <- list( c("1 Bacillus grandii", "2 Bacillus atticus"), c("3 Bacillus rossius", "2 Bacillus atticus"), c("1 Bacillus grandii", "3 Bacillus rossius") )
m_module_dNdS_bgm <- subset(male_gonad_network_tot, BGM_dN<3 & BGM_dNdS<1)
m_module_dNdS_bgm <- m_module_dNdS_bgm$BGM_dNdS
m_module_dNdS_bgm <- data.frame(group = "1 Bacillus grandii", value = m_module_dNdS_bgm)
m_module_dNdS_bat <- subset(male_gonad_network_tot, BAT_dN<3 & BAT_dNdS<1)
m_module_dNdS_bat <- m_module_dNdS_bat$BAT_dNdS
m_module_dNdS_bat <- data.frame(group = "2 Bacillus atticus", value = m_module_dNdS_bat)
m_module_dNdS_bro <- subset(male_gonad_network_tot, BRO_dN<3 & BRO_dNdS<1)
m_module_dNdS_bro <- m_module_dNdS_bro$BRO_dNdS
m_module_dNdS_bro <- data.frame(group = "3 Bacillus rossius", value = m_module_dNdS_bro)
m_module_dNdS_plot.data <- rbind(m_module_dNdS_bgm,m_module_dNdS_bat,m_module_dNdS_bro)
m_module_dNdS_plot <- ggplot(m_module_dNdS_plot.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
labs(x=" ", y = "dNdS", title = "genes in modules associated with B.grandii male gonad") + theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
scale_fill_manual(values=c("#04b4d6", "#ffc037", "#f55a4f")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
f_module_dNdS_bgm <- subset(female_gonad_network_tot, BGM_dN<3 & BGM_dNdS<1)
f_module_dNdS_bgm <- f_module_dNdS_bgm$BGM_dNdS
f_module_dNdS_bgm <- data.frame(group = "1 Bacillus grandii", value = f_module_dNdS_bgm)
f_module_dNdS_bat <- subset(female_gonad_network_tot, BAT_dN<3 & BAT_dNdS<1)
f_module_dNdS_bat <- f_module_dNdS_bat$BAT_dNdS
f_module_dNdS_bat <- data.frame(group = "2 Bacillus atticus", value = f_module_dNdS_bat)
f_module_dNdS_bro <- subset(female_gonad_network_tot, BRO_dN<3 & BRO_dNdS<1)
f_module_dNdS_bro <- f_module_dNdS_bro$BRO_dNdS
f_module_dNdS_bro <- data.frame(group = "3 Bacillus rossius", value = f_module_dNdS_bro)
f_module_dNdS_plot.data <- rbind(f_module_dNdS_bgm,f_module_dNdS_bat,f_module_dNdS_bro)
f_module_dNdS_plot <- ggplot(f_module_dNdS_plot.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
labs(x=" ", y = "dNdS", title = "genes in modules associated with B.grandii female gonad") + theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
scale_fill_manual(values=c("#04b4d6", "#ffc037", "#f55a4f")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
nonmodule_dNdS_bgm <- subset(all_other_network_tot, BGM_dN<3 & BGM_dNdS<1)
nonmodule_dNdS_bgm <- nonmodule_dNdS_bgm$BGM_dNdS
nonmodule_dNdS_bgm <- data.frame(group = "1 Bacillus grandii", value = nonmodule_dNdS_bgm)
nonmodule_dNdS_bat <- subset(all_other_network_tot, BAT_dN<3 & BAT_dNdS<1)
nonmodule_dNdS_bat <- nonmodule_dNdS_bat$BAT_dNdS
nonmodule_dNdS_bat <- data.frame(group = "2 Bacillus atticus", value = nonmodule_dNdS_bat)
nonmodule_dNdS_bro <- subset(all_other_network_tot, BRO_dN<3 & BRO_dNdS<1)
nonmodule_dNdS_bro <- nonmodule_dNdS_bro$BRO_dNdS
nonmodule_dNdS_bro <- data.frame(group = "3 Bacillus rossius", value = nonmodule_dNdS_bro)
nonmodule_dNdS_plot.data <- rbind(nonmodule_dNdS_bgm,nonmodule_dNdS_bat,nonmodule_dNdS_bro)
nonmodule_dNdS_plot <- ggplot(nonmodule_dNdS_plot.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
labs(x=" ", y = "dNdS", title = "genes non related to B. grandii gonad modules")+ theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
scale_fill_manual(values=c("#04b4d6", "#ffc037", "#f55a4f")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
###
my_comparisons <- list( c("female gonads", "male gonads"), c("gonads unrelated", "male gonads"), c("gonads unrelated", "female gonads") )
#
(shapiro.test(m_module_dNdS_bgm$value))$p.value## [1] 3.51044e-33(shapiro.test(f_module_dNdS_bgm$value))$p.value## [1] 4.716072e-23(shapiro.test(nonmodule_dNdS_bgm$value))$p.value## [1] 1.973568e-40(shapiro.test(m_module_dNdS_bat$value))$p.value## [1] 6.134926e-34(shapiro.test(f_module_dNdS_bat$value))$p.value## [1] 2.097348e-23(shapiro.test(nonmodule_dNdS_bat$value))$p.value## [1] 9.356442e-40(shapiro.test(nonmodule_dNdS_bro$value))$p.value## [1] 1.278275e-41(shapiro.test(nonmodule_dNdS_bro$value))$p.value## [1] 1.278275e-41(shapiro.test(nonmodule_dNdS_bro$value))$p.value## [1] 1.278275e-41m_module_dNdS_bgm <- subset(male_gonad_network_tot, BGM_dN<3 & BGM_dNdS<1)
m_module_dNdS_bgm <- m_module_dNdS_bgm$BGM_dNdS
m_module_dNdS_bgm <- data.frame(group = "m", value = m_module_dNdS_bgm)
m_module_dNdS_bat <- subset(male_gonad_network_tot, BAT_dS>0.001 & BAT_dN<3 & BAT_dNdS<1)
m_module_dNdS_bat <- m_module_dNdS_bat$BAT_dNdS
m_module_dNdS_bat <- data.frame(group = "m", value = m_module_dNdS_bat)
m_module_dNdS_bro <- subset(male_gonad_network_tot,BRO_dN<3 & BRO_dNdS<1)
m_module_dNdS_bro <- m_module_dNdS_bro$BRO_dNdS
m_module_dNdS_bro <- data.frame(group = "m", value = m_module_dNdS_bro)
f_module_dNdS_bgm <- subset(female_gonad_network_tot, BGM_dN<3 & BGM_dNdS<1)
f_module_dNdS_bgm <- f_module_dNdS_bgm$BGM_dNdS
f_module_dNdS_bgm <- data.frame(group = "f", value = f_module_dNdS_bgm)
f_module_dNdS_bat <- subset(female_gonad_network_tot,BAT_dN<3 & BAT_dNdS<1)
f_module_dNdS_bat <- f_module_dNdS_bat$BAT_dNdS
f_module_dNdS_bat <- data.frame(group = "f", value = f_module_dNdS_bat)
f_module_dNdS_bro <- subset(female_gonad_network_tot, BRO_dN<3 & BRO_dNdS<1)
f_module_dNdS_bro <- f_module_dNdS_bro$BRO_dNdS
f_module_dNdS_bro <- data.frame(group = "f", value = f_module_dNdS_bro)
nonmodule_dNdS_bgm <- subset(all_other_network_tot, BGM_dN<3 & BGM_dNdS<1)
nonmodule_dNdS_bgm <- nonmodule_dNdS_bgm$BGM_dNdS
nonmodule_dNdS_bgm <- data.frame(group = "o", value = nonmodule_dNdS_bgm)
nonmodule_dNdS_bat <- subset(all_other_network_tot, BAT_dN<3 & BAT_dNdS<1)
nonmodule_dNdS_bat <- nonmodule_dNdS_bat$BAT_dNdS
nonmodule_dNdS_bat <- data.frame(group = "o", value = nonmodule_dNdS_bat)
nonmodule_dNdS_bro <- subset(all_other_network_tot, BRO_dN<3 & BRO_dNdS<1)
nonmodule_dNdS_bro <- nonmodule_dNdS_bro$BRO_dNdS
nonmodule_dNdS_bro <- data.frame(group = "o", value = nonmodule_dNdS_bro)
dNdS_plot_bro.data <- rbind(m_module_dNdS_bro,f_module_dNdS_bro,nonmodule_dNdS_bro)
dNdS_plot_bro.data[dNdS_plot_bro.data == "m"] <- "male gonads"
dNdS_plot_bro.data[dNdS_plot_bro.data == "f"] <- "female gonads"
dNdS_plot_bro.data[dNdS_plot_bro.data == "o"] <- "gonads unrelated"
dNdS_plot_bro <- ggplot(dNdS_plot_bro.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
labs(x=" ", y = "dNdS", title = "Bacillus rossius")+ theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
scale_fill_manual(values=c("#f55a4f", "#f55a4f", "#f55a4f")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
dNdS_plot_bat.data <- rbind(m_module_dNdS_bat,f_module_dNdS_bat,nonmodule_dNdS_bat)
dNdS_plot_bat.data[dNdS_plot_bat.data == "m"] <- "male gonads"
dNdS_plot_bat.data[dNdS_plot_bat.data == "f"] <- "female gonads"
dNdS_plot_bat.data[dNdS_plot_bat.data == "o"] <- "gonads unrelated"
dNdS_plot_bat <- ggplot(dNdS_plot_bat.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
labs(x=" ", y = "dNdS", title = "Bacillus atticus")+ theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
scale_fill_manual(values=c("#ffc037", "#ffc037", "#ffc037")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
dNdS_plot_bgm.data <- rbind(m_module_dNdS_bgm,f_module_dNdS_bgm,nonmodule_dNdS_bgm)
dNdS_plot_bgm.data[dNdS_plot_bgm.data == "m"] <- "male gonads"
dNdS_plot_bgm.data[dNdS_plot_bgm.data == "f"] <- "female gonads"
dNdS_plot_bgm.data[dNdS_plot_bgm.data == "o"] <- "gonads unrelated"
dNdS_plot_bgm <- ggplot(dNdS_plot_bgm.data, aes(x=group, y=value, fill=group)) + geom_boxplot(width=0.2, outlier.size=0, alpha = 0.7, outlier.shape=NA) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif", label.y = c(0.6,0.7,0.8)) +
labs(x=" ", y = "dNdS", title = "Bacillus grandii")+ theme_classic() +
theme(axis.text=element_text(size=10), axis.title=element_text(size=12,face="bold"), title=element_text(size=16,face="bold")) + ylim(0, 1) +
scale_fill_manual(values=c("#04b4d6", "#04b4d6", "#04b4d6")) + theme(legend.position = "none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
coord_fixed(ratio=3)
###
grid.arrange(nonmodule_dNdS_plot,
f_module_dNdS_plot,
m_module_dNdS_plot,
dNdS_plot_bro,
dNdS_plot_bat,
dNdS_plot_bgm,
ncol=3)
calculate B. atticus gene connectivity
adj_bgm <- adjacency(datExpr0, selectCols = NULL, type = "signed",
power = 19,
corFnc = "cor", corOptions = "use = 'p', method = 'spearman'",
distFnc = "dist", distOptions = "method = 'euclidean'")
intramodularConnectivity_bgm <- intramodularConnectivity(adj_bgm, mergedColors, scaleByMax = FALSE)
intramodularConnectivity_bgm['OG'] <- row.names(intramodularConnectivity_bgm)
data wrangling pt.2
tot <- read.table(file = "def.tab", sep = " ", header= TRUE)
yellow_tot <- subset(tot, tot$OG %in% yellow_module)
green_tot <- subset(tot, tot$OG %in% green_module)
turquoise_tot <- subset(tot, tot$OG %in% turquoise_module)
red_tot <- subset(tot, tot$OG %in% red_module)
brown_tot <- subset(tot, tot$OG %in% brown_module)
blue_tot <- subset(tot, tot$OG %in% blue_module)
black_tot <- subset(tot, tot$OG %in% black_module)
#
male_gonad_network_tot <- rbind(yellow_tot,green_tot,turquoise_tot)
male_gonad_network_genes <- male_gonad_network_tot$OG
female_gonad_network_tot <- rbind(red_tot,brown_tot)
female_gonad_network_genes <- female_gonad_network_tot$OG
all_other_network_tot <- subset(tot, !(tot$OG %in% male_gonad_network_genes))
all_other_network_tot <- subset(all_other_network_tot, !(tot$OG %in% female_gonad_network_genes))
all_other_network_genes <- all_other_network_tot$OG
#
con_tot <- merge(intramodularConnectivity_bgm, tot, by="OG")
names(con_tot)[names(con_tot) == "kTotal"] <- "kTotal_bgm"
names(con_tot)[names(con_tot) == "kWithin"] <- "kWithin_bgm"
names(con_tot)[names(con_tot) == "kOut"] <- "kOut_bgm"
names(con_tot)[names(con_tot) == "kDiff"] <- "kDiff_bgm"
con_tot_mod_other <- subset(con_tot, con_tot$OG %in% all_other_network_genes)
con_tot_DE <- subset(con_tot, con_tot$OG %in% male_gonad_DE_genes)
con_tot_DE_other <- subset(con_tot, !(con_tot$OG %in% male_gonad_DE_genes))
mg_module_genes <- c(turquoise_module,green_module,yellow_module)
mg_module_genes_con_tot <- subset(con_tot, tot$OG %in% mg_module_genes)
fg_module_genes <- c(brown_module, red_module)
fg_module_genes_con_tot <- subset(con_tot, tot$OG %in% fg_module_genes)
gene connectivity VS LogFC
# bgm VS bat
bgm_kTotal_bat_LogFC_spearman_mg_modules <- cor.test(mg_module_genes_con_tot$kTotal_bgm , mg_module_genes_con_tot$BAT_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bat_LogFC_plot_mg_modules <- ggplot() +
geom_point(data=mg_module_genes_con_tot, aes(kTotal_bgm, BAT_logFC, color=BAT_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. atticus (asex)", title = paste0(
"genes in male gonad modules pval: ", round(bgm_kTotal_bat_LogFC_spearman_mg_modules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bat_LogFC_spearman_mg_modules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#ffc037", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bat_LogFC_spearman_fg_modules <- cor.test(fg_module_genes_con_tot$kTotal_bgm , fg_module_genes_con_tot$BAT_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bat_LogFC_plot_fg_modules <- ggplot() +
geom_point(data=fg_module_genes_con_tot, aes(kTotal_bgm, BAT_logFC, color=BAT_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. atticus (asex)", title = paste0(
"genes in female gonad modules pval: ", round(bgm_kTotal_bat_LogFC_spearman_fg_modules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bat_LogFC_spearman_fg_modules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#ffc037", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bat_LogFC_spearman_nonmodules <- cor.test(con_tot_mod_other$kTotal_bgm , con_tot_mod_other$BAT_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bat_LogFC_plot_nonmodules <- ggplot() +
geom_point(data=con_tot_mod_other, aes(kTotal_bgm, BAT_logFC,color=BAT_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. atticus (asex)", title = paste0(
"all other genes pval: ", round(bgm_kTotal_bat_LogFC_spearman_nonmodules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bat_LogFC_spearman_nonmodules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#ffc037", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
# bgm VS bro
bgm_kTotal_bro_LogFC_spearman_mg_modules <- cor.test(mg_module_genes_con_tot$kTotal_bgm , mg_module_genes_con_tot$BRO_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bro_LogFC_plot_mg_modules <- ggplot() +
geom_point(data=mg_module_genes_con_tot, aes(kTotal_bgm, BRO_logFC, color=BRO_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. rossius (asex)", title = paste0(
"genes in male gonad modules pval: ", round(bgm_kTotal_bro_LogFC_spearman_mg_modules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bro_LogFC_spearman_mg_modules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#f55a4f", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bro_LogFC_spearman_fg_modules <- cor.test(fg_module_genes_con_tot$kTotal_bgm , fg_module_genes_con_tot$BRO_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bro_LogFC_plot_fg_modules <- ggplot() +
geom_point(data=fg_module_genes_con_tot, aes(kTotal_bgm, BRO_logFC, color=BRO_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. rossius (asex)", title = paste0(
"genes in female gonad modules pval: ", round(bgm_kTotal_bro_LogFC_spearman_fg_modules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bro_LogFC_spearman_fg_modules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#f55a4f", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bro_LogFC_spearman_nonmodules <- cor.test(con_tot_mod_other$kTotal_bgm , con_tot_mod_other$BRO_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bro_LogFC_plot_nonmodules <- ggplot() +
geom_point(data=con_tot_mod_other, aes(kTotal_bgm, BRO_logFC, color=BRO_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. rossius (asex)", title = paste0(
"all other genes pval: ", round(bgm_kTotal_bro_LogFC_spearman_nonmodules$p.value,digits = 3), " rho: ", round(bgm_kTotal_bro_LogFC_spearman_nonmodules$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#f55a4f","white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
# bgm VS bgm F
bgm_kTotal_bgm_LogFC_spearman_mg_modules_f <- cor.test(mg_module_genes_con_tot$kTotal_bgm , mg_module_genes_con_tot$BGM_F_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_mg_modules_f <- ggplot() +
geom_point(data=mg_module_genes_con_tot, aes(kTotal_bgm, BGM_F_logFC, color=BGM_F_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii F (sex)", title = paste0(
"genes in male gonad modules pval: ", round(bgm_kTotal_bgm_LogFC_spearman_mg_modules_f$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_mg_modules_f$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bgm_LogFC_spearman_fg_modules_f <- cor.test(fg_module_genes_con_tot$kTotal_bgm , fg_module_genes_con_tot$BGM_F_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_fg_modules_f <- ggplot() +
geom_point(data=fg_module_genes_con_tot, aes(kTotal_bgm, BGM_F_logFC, color=BGM_F_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii F (sex)", title = paste0(
"genes in female gonad modules pval: ", round(bgm_kTotal_bgm_LogFC_spearman_fg_modules_f$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_fg_modules_f$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bgm_LogFC_spearman_nonmodules_f <- cor.test(con_tot_mod_other$kTotal_bgm , con_tot_mod_other$BGM_F_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_nonmodules_f <- ggplot() +
geom_point(data=con_tot_mod_other, aes(kTotal_bgm, BGM_F_logFC, color=BGM_F_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii F (sex)", title = paste0(
"all other genes pval: ", round(bgm_kTotal_bgm_LogFC_spearman_nonmodules_f$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_nonmodules_f$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
# bgm VS bgm M
bgm_kTotal_bgm_LogFC_spearman_mg_modules_m <- cor.test(mg_module_genes_con_tot$kTotal_bgm , mg_module_genes_con_tot$BGM_M_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_mg_modules_m <- ggplot() +
geom_point(data=mg_module_genes_con_tot, aes(kTotal_bgm, BGM_M_logFC, color=BGM_M_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii M (sex)", title = paste0(
"genes in male gonad modules pval: ", round(bgm_kTotal_bgm_LogFC_spearman_mg_modules_m$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_mg_modules_m$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bgm_LogFC_spearman_fg_modules_m <- cor.test(fg_module_genes_con_tot$kTotal_bgm , fg_module_genes_con_tot$BGM_M_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_fg_modules_m <- ggplot() +
geom_point(data=fg_module_genes_con_tot, aes(kTotal_bgm, BGM_M_logFC, color=BGM_M_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii M (sex)", title = paste0(
"genes in female gonad modules pval: ", round(bgm_kTotal_bgm_LogFC_spearman_fg_modules_m$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_fg_modules_m$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
bgm_kTotal_bgm_LogFC_spearman_nonmodules_m <- cor.test(con_tot_mod_other$kTotal_bgm , con_tot_mod_other$BGM_M_logFC, method = "spearman", continuity = FALSE, conf.level = 0.95)
bgm_kTotal_bgm_LogFC_plot_nonmodules_m <- ggplot() +
geom_point(data=con_tot_mod_other, aes(kTotal_bgm, BGM_M_logFC, color=BGM_M_padj), size=5, alpha = 0.5, shape= 20) +
labs(x="total connectivity B. grandii (sex)", y = "LogFC B. grandii M (sex)", title = paste0(
"all other genes pval: ", round(bgm_kTotal_bgm_LogFC_spearman_nonmodules_m$p.value,digits = 3), " rho: ", round(bgm_kTotal_bgm_LogFC_spearman_nonmodules_m$estimate,digits = 3))) +
theme_classic() + theme(axis.text=element_text(size=12), axis.title=element_text(size=16,face="bold")) +
geom_hline(yintercept=0, color = "black",size=0.5, alpha = 0.2) +
ylim(-7, +7) + scale_colour_gradientn(colours = c("#04b4d6", "white","white","white","white"), breaks = 0.01, name = "padj") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(colour = "black", size=1)) +
theme(plot.margin = unit(c(1,1,1,1), "cm")) + theme(plot.title = element_text(size=16,face="bold"))
grid.arrange(bgm_kTotal_bgm_LogFC_plot_mg_modules_f,
bgm_kTotal_bgm_LogFC_plot_mg_modules_m,
bgm_kTotal_bat_LogFC_plot_mg_modules,
bgm_kTotal_bro_LogFC_plot_mg_modules,
bgm_kTotal_bgm_LogFC_plot_fg_modules_f,
bgm_kTotal_bgm_LogFC_plot_fg_modules_m,
bgm_kTotal_bat_LogFC_plot_fg_modules,
bgm_kTotal_bro_LogFC_plot_fg_modules,
bgm_kTotal_bgm_LogFC_plot_nonmodules_f,
bgm_kTotal_bgm_LogFC_plot_nonmodules_m,
bgm_kTotal_bat_LogFC_plot_nonmodules,
bgm_kTotal_bro_LogFC_plot_nonmodules,
layout_matrix = rbind(c(1,2,3,4),c(5,6,7,8),c(9,10,11,12)))